Testing Skyrme energy-density functionals with the QRPA in low-lying vibrational 

states of rare-earth nuclei 
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Although nuclear energy density functionals are determined primarily by fitting to ground state 
properties, they are often applied in nuclear astrophysics to excited states, usually through the 
quasiparticle random phase approximation (QRPA). Here we test the Skyrme functionals SkM* and 
SLy4 along with the self-consistent QRPA by calculating properties of low-lying vibrational states 
in a large number of well-deformed even-even rare-earth nuclei. We reproduce trends in energies 
and transition probabilities associated with 7- vibrational states, but our results are not perfect and 
indicate the presences of multi-particle-hole correlations that are not included in the QRPA. The 
Skyrme functional SkM* performs noticeably better than SLy4. In a few nuclei, changes in the 
treatment of the pairing energy functional have a significant effect. The QRPA is less successful 
with "/^-vibrational" states than with the 7-vibrational states. 

PACS numbers: 21.10.Re, 21.60.Jz, 27.70.-hq 



I. INTRODUCTION 

Modern supercomputers are making the quantitative 
theoretical treatment of nuclear structure increasingly 
common. In light nuclei, Greens-function Monte Carlo 
methods [1, 2] and the no-core shell model [3, 4] yield ac- 
curate ab initio results, and in medium-mass nuclei the 
coupled cluster method [5, 6] is proving successful. In nu- 
clei with A > 50, techniques related to density-functional 
theory (DFT) [7] are the state of the art. Accuracy, at 
least for ground-state properties, is limited only by the 
quality of the functionals, which are continually improv- 
ing [8]. 

One advantage of DFT is its applicability to nearly all 
heavy nuclei. Such flexibility is particularly important 
for nuclear astrophysics, which attempts to explain the 
synthesis of all the elements. Another advantage is a 
natural extension, through the self-consistent quasiparti- 
cle random phase approximation (QRPA) to excitations. 
Excited states are as important as ground states in many 
nucleosynthetic reactions, and so Goriely et al. [9], for 
example, used the QRPA to compute radiative neutron 
capture in a wide range of nuclei. Such calculations, how- 
ever, have generally ignored deformation, or treated it in 
a crude way. The logical next step is to take the effects 
of deformation into account in a self-consistent fashion. 

Fortunately, self-consistent QRPA calculations in 
heavy deformed nuclei are now becoming possible. Re- 
cently, we developed a scaled parallel Skyrme-QRPA 
code [10] for arbitrary axially-deformed (parity conserv- 
ing) even-even nuclei. Our code is one of the few [11] 
to treat heavy deformed nuclei in the QRPA without 
simplification. (For other calculations, including those 
in lighter nuclei and those in the RPA, the spherical 
QRPA, and separable approximations, see the work cited 



in Ref. [10] and, e.g., the more recent Ref. [12].) In 
this paper, we present calculations with two Skyrme 
energy-density functionals of properties of low-energy vi- 
brational states in rare-earth nuclei. As promised in Ref. 
[10], we discuss the performance of both the functionals 
and the QRPA. 

In Sec. II below we list the nuclei that we explore 
and present technical information about our calculations. 
In Sec. Ill we show results for 7-vibrational states and 
discuss the performance of the Skyrme QRPA, which 
we compare with methods used in earlier calculations. 
Sec IV treats "^^-vibrational" states-^ briefly, and Sec. V 
is a conclusion. An appendix presents equations for two- 
body matrix elements of the Coulomb- direct interaction 
and discusses computational efliciency. 



II. SELECTION OF NUCLEI AND METHOD 
OF CALCULATION 

The vibrational states we examine are all in rare-earth 
nuclei. The advantage of this region of the isotopic chart 
is the abundance of reliable experimental data [14-42], 
accumulated over the last half century. Multiple results 
exist for many of the nuclei and there are few serious dis- 
crepancies. In addition, the large deformation of many 
of the rare earths make them better candidates for a suc- 
cessful QRPA treatment than transitional nuclei, which 
tend to be soft. We choose the 27 nuclei shown in Fig. 1 
for our calculation. They are all axially symmetric and 
well deformed, with f3 > 0.3, and for all but a few the 
energies of their 7- vibrations {K^ = 2+, the second or 
third = 2+ states) have been measured and appear 
in Ref. [43]. We calculate 7-vibrational energies and E2 
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^ We put the term in quotes to indicate that many of those states 
are not purely vibrational [13]. 
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FIG. 1: (Color online) Rare-earth region of the isotopic chart. 
Shaded area shows nuclei with deformation /3 > 0.3 in HFB 
calculations with SkM* (unpublished, see Ref. [47] for SLy4 
which gives a similar result). We calculate energies and tran- 
sition probabilities of 7-vibrational states for all nuclei whose 
isotopic symbols appear in the figure. Squares without sym- 
bols correspond to nuclei for which experimental data on 7 
vibrations are not in Ref. [43]. 



suit. In the other rare-earth nuclei, the dimension ranges 
from 19000 to 28000. 

Spurious states associated with particle-number con- 
servation make the necessary space much larger for 
"/3-vibrations." There we use ((^pa^r)^ ((^ph*)^) = 
(10-^^,10-^) with Ecut = 200 MeV, values with which 
the dimension of the two-canonical-quasiparticle space is 
60000 to 75000. Even with this large dimension, how- 
ever, the spurious state does not separate perfectly and 
we present results only for cases in which the separation 
is good. 

Deriving the QRPA equations for an axially-symmetric 
system is tedious but not difficult and can be done by 
starting from the general equations in, e.g., Ref. [46]. 
In the appendix, therefore, we display only our repre- 
sentation of the Coulomb-direct matrix elements. These 
require more numerical effort than matrix elements of a 
(^-interaction, and so benefit more from a computation- 
ally efficient procedure. 



III. 7-VIBRATIONS 



excitation strengths in all 27 nuclei with the Skyrme func- 
tionals SkM* [44] and SLy4 [45], and in a few nuclei 
we do the same for vibrational" states {K^ = 0+) 
with SkM*. We use the traditional volume-pairing en- 
ergy functional [46] for simplicity. 

Our procedure has two steps: a Hartree-Fock- 
Bogoliubov (HFB) calculation with the Vanderbilt HFB 
code [48], and a QRPA calculation that uses the results 
of the HFB run. Both steps use B splines [49-51] to rep- 
resent wave functions on a 42 by 42 cylindrical mesh with 
< z^p < 20 fm. We use box boundary conditions to 
discretize the continuum, and introduce a quasiparticle 
cutoff energy E^ut of 60 MeV or 200 MeV in the HFB cal- 
culation to limit the set of quasiparticle wave functions 
that determine the density and the pairing tensor. The 
two cutoffs require different pairing strengths, which we 
adjust via the three-point formula [52] so as to reproduce 
the pairing gaps of ^^^Yb (obtained from experimental 
masses). In the other nuclei, this procedure usually re- 
produces overall pairing gaps to within ±150 keV. We 
restrict the z-component of the angular momentum of 
the wave functions to be less than or equal to 19/2 h. 

Next we transform the quasiparticle wave functions 
to the canonical-basis and introduce two cutoff occupa- 
tion probabilities ('^pai^)^ ('^ph*)^' used also in our 
prior work [10, 46, 53], to truncate the two-canonical- 
quasiparticle basis in which we construct the QRPA 
Hamiltonian matrix. We take (('^pair)^? (('^ph*)^) = 
(10-^^,10-^) for Ec^t = 60 MeV, and (10-^ 10"^) for 
^cut = 200 MeV in the 7-vibration calculation with 
SkM*. Those values make the dimension of the two- 
canonical-quasiparticle basis about 22000 in ^^^Yb, a 
number that is large enough to yield a convergent re- 



A. Energies and transition strengths 

Figure 2 shows measured 7-vibration energies along- 
side the results of our two QRPA calculations and 
a collective-model calculation (with parameters deter- 
mined from the Gogny energy functional) by Delaroche 
et al. [54]. In all the plots the minimum energy occurs 
around ^4=162. The minimum in the Dy and Er isotopes 
is at A/" = 98, both in the data and the SkM*. The Yb 
isotopes are particularly well reproduced by the SkM* 
calculation. But the QRPA calculations show a stronger 
A-dependence than the data, with the SLy4 results show- 
ing the strongest dependence. And overall, neither of the 
QRPA calculations is as good as that of Ref. [54] . 

Figure 3 shows E2 transition strengths B{E2; 0+ 
2+), hereafter denoted B{E2)'\', for the same isotopes. 
Overall, the calculations reproduce the data reasonably 
well, except in ^^^Dy with SLy4, and again are partic- 
ularly good in the Yb isotopes. As before, SkM* is no- 
ticeably better than SLy4. The energies and B{E2)Ys in 
our calculations are anticorrelated in general, a feature 
expected of harmonic vibrations. On the other hand, the 
experimental 5(£^2)t's in Er decrease monotonically with 
A, even though the dependence of the energy is slightly 
parabolic. 

To characterize the performance of the two functionals 
statistically, we introduce, following Refs. [53, 55] the 
measures 

=ln(^cal/^exp) (1) 

and 

Rq = In ^5(^2)tcal /5(^2)texp , (2) 
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FIG. 3: (Color online) B(^2;0+ 2+) corresponding to Fig. 2. The value for ^^^Dy in c) is 0.562 e^b^ This figure has no 
panel d) because the results from the calculation of Delaroche et al. [54] are not published. We include only those experimental 
data that are labeled 7- vibrations in Ref . [43] . The symbols for particular isotopic chains are the same in each panel. 



where Ecai and £^exp are the calculated and experimen- 
tal energies of the 7-vibrational state. The results are 
in Tab. L SLy4 actually does better than SkM* in the 
averages, but gives much larger dispersions. 

Table II shows the statistical measures for the spherical 



nuclei treated in Ref. [53] and for the subset of those 
nuclei that exhibit "low softness." (Some of the other 
nuclei in Ref. [53] are transitional.) There are far more 
nuclei in the spherical data set than in the deformed rare- 
earth set, so it is hard to make a precise comparison of 
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TABLE I: Average of Re (Re), dispersion of Re {cte), and 
the same for Rq. 





Re 


aE 


Rq 




SkM* 


0.28 


0.18 


-0.13 


0.14 


SLy4 


0.20 


0.50 


-0.004 


0.31 



TABLE IL Statistical measures for the spherical nuclei and 
for the subset with low softness from Ref. [53]. Rq and aq 
were not calculated separately for the low-softness nuclei in 
that paper because of a lack of E2 data. 







Re 


(TE 


Rq 




All 


SkM* 


0.11 


0.44 


-0.29 


0.53 


SLy4 


0.33 


0.51 


-0.32 


0.42 


Low Softness 


SkM* 


0.27 


0.35 








SLy4 


0.47 


0.48 







performance. But deformation does not appear to affect 
it significantly. 



B. N- and Z-dependence 

Our calculations show a stronger dependence on N 
than do the data in most isotopic chains, behavior that 
may be due to insufficient configuration mixing in our cal- 
culations. Figure 4 shows the TV-dependence of the calcu- 
lated 7-vibrational energy, the two-quasiparticle energy 
-^2qp of component, in the quasiparticle basis, with 
the largest neutron forward amplitude, and the absolute 
value of the backward amplitude |lxn| of the same com- 
ponent, all for SkM*. (We transformed amplitudes from 
the canonical-quasiparticle basis to do this analysis.) For 
N < 100, The 7-vibrational energy is positively corre- 
lated with £^2^p, and anticorrelated with llxnl, indicating 
a connection between the A^-dependence of those solu- 
tions and a particular two-quasiparticle state. The down- 
ward shift of about 1 MeV between the two-quasiparticle 
energy and the full QRPA energy, seen in panels a) and 
b), then characterizes the effect of the residual interac- 
tion. Fig. 5 shows all the same phenomena in the Z 
dependence of our results, except in the energies of the 
N = 102 isotones. 

Figure 6 shows the absolute values of the nine largest 
neutron forward amplitudes in three Dy isotopes around 
i^^Dy, which is the one with the lowest phonon energy. 
Clearly the two largest components are far more impor- 
tant than the rest. And though we don't show it, a 
similar curve characterizes the protons. From all this, 
we conclude that the two-quasiparticle state with the 
largest neutron forward amplitudes plays a significant 
role in the A/'-dependences of the QRPA solutions, and 
that the same statement is true of proton forward am- 
plitudes and Z dependence. (The second largest com- 



ponents are potentially also important.) The weaker N- 
and Z-dependence in the data suggests that we exagger- 
ate the importance of those particular two-quasiparticle 
states, perhaps by underestimating configuration mixing. 
It is quite possible that a better solution requires many- 
body correlations beyond the QRPA. 

Figure 7 shows ^2qp SLy4 calculation. Inter- 

estingly, the range of the £'2qp close to that produced 
by SkM*, as one can see by comparing with panel b) 
of Fig. 4. We conclude that the effects of the residual 
interaction on A dependence are quite different in the 
two calculations, leading to the noticeable differences in 
Fig. 2. 



C. ^-pairing functional 

Low-energy quasiparticle states are obviously affected 
by the choice of pairing functional, and the volume pair- 
ing we use can be varied without worsening its ability to 
reproduce pairing gaps. The reason is that the ^-function 
interaction is singular in the pairing channel and so must 
be regularized (see, e. g. Ref. [56]). Here we do so by 
cutting off the single-quasiparticle spectrum. This pro- 
cedure makes the strength of the interaction depend on 
the cutoff as well as on the experimental pairing gaps to 
which it is fit. To illustrate the effect of the cutoff on 7 
vibrations, we show in Fig. 8 the results of calculations 
with the two different cutoffs £^cut mentioned in Sec. II. 
The two cutoffs require different pairing strengths Gq (for 
the values, see the caption) to ensure similar predictions 
for pairing gaps. We refer to the two calculations as A 
and B, with A having the smaller E^cut, and therefore the 
larger pairing strength. 

The differences in the results are mostly minor, but the 
5(E2)t in i^^Dy, the energy and B(£'2)t in ^^^Hf, and 
the energy in ^''^Yb are all quite different. We account 
for the differences in ^^^Dy by referring to Fig. 9. Since 
calculation A uses a larger pairing strength, it produces 
higher energies for low-lying quasiparticles than does cal- 
culation B. In the separable approximation [57, 58], the 
forward QRPA amplitudes can be written as 



^ qpl,qp2 OC 1/ (E'qpl + £'qp2 " ^7) , 



(3) 



where qpl and qp2 denote quasiparticle states, and Ej 
is the energy of the 7-vibrational state. Using Eq. (3) 
and the values read from the figures, one can estimate 
the ratio of the forward amplitudes of the largest two- 
quasiprarticle component in calculations A and B. The 
result, under the assumption that the interaction matrix 
elements are the same in the two calculations, is 



^^qpl,qp2 



^qpl,qp2 



C± 0.9. 



(4) 



Panel b of Fig. 9 shows that the exact ratio is 0.8, so that 
half the difference between the two calculations can be 
explained by considering only the quasiparticle energies. 
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FIG. 5: (Color online) The same panels as Fig. 4 but for protons rather than neutrons. The connected points are isotones, and 
the same symbol indicates a given isotone in each panel, and |lxp| now refer to the proton components. 



This analysis implies that the QRPA solution is sensitive 
to the energies of important quasiparticle states, and thus 
to the pairing functional, when those energies are small. 
One can take advantage of this to fix the cutoff energy 
as well as the pairing strength by fitting to properties 
that depend sensitvely on low-energy quasiparticle states. 
Our calculation shows, for example, that Ecut = 60 MeV 
is better than 200 MeV. 

The differences between calculations A and B in ^^^Yb 



and ^^^Hf are more complicated (the corresponding two- 
quasiparticle components do not have the same order as 
Fig. 9), and we could not find a simple explanation for 
them. And changes in the pairing cutoff are clearly not 
enough to fix the problem with A^- and Z-dependence in 
the QRPA. 
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FIG. 6: (Color online) Absolute values of the nine largest 
neutron forward amplitudes in 162-166-^^ rpj^^ integer on the 
X-axis labels the two-quasiparticle components. 



D. Comparison with older calculations of 
7-vibrational states 

Early work on vibrations in rare-earth nuclei of- 
ten made use of the pairing-plus-QQ (quadrupole- 
quadrupole) Hamiltonian, both in the (Q)RPA [59-63] 
and in approximations that went beyond the QRPA or- 
der, e. g. [64-66]. Single-particle energies were usually 
obtained from the Nilsson potential, with slight shifts to 
improve phenomenology, and the strength of the QQ in- 
teractions was modified slightly from the self-consistent 
value so as to reproduce the energies of the 7-vibrational 
states. The adjustment to energies means that B{E2)']''s 
are an important test of the model's predictive power. 

Figure 10 shows the energies and B{E2) ^ from 
Ref. [60]. Their energies were perhaps not quite as good 
overall as ours, but also did not exhibit the sharp min- 
imum we get around A ^ 164. The authors themselves 
stated that no single interaction strength reproduces the 
energies of all nuclei calculated. Their B{E2)"\'^s are too 
large by a factor of two or more, a deficiency that was 
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FIG. 7: (Color online) ^^^p produced by SLy4. The symbols 
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pointed out again in Refs. [59, 62]. Marshalek et al. [59] 
listed approximations that might cause problems in pre- 
dicted B(^2)t's. Since we do better in that observable, 
we believe that the cause is in fact the interaction. 

Rare-earth 7 vibrations have also been addressed in 
other models. Reference [64], using the boson-expansion 
method, obtained B{E2)^ = 0.130 e^b^ in ^^^^Sm, a value 
somewhat larger than ours. Soloviev et al. [65, 66] used 
the quasiparticle-phonon nuclear model, which includes 
two-phonon couplings, and obtained B{E2) ^ = 0.127 
g2^2 (168^^)^ Q g2^2 (iT2Yb), and 0.122 e^b^ {^'^m). 

Those transition probabilities are close to the experimen- 
tal data (0.116 e^b^in ^^^Hf [65]), and the fit of the in- 
teraction meant that energies were also reproduced well. 
See also Ref. [63] which used a modified QQ and an ef- 
fective three-body interactions. 

Explicitly collective models have also been used. Ku- 
mar [67] obtained an energy of 1.438 MeV (close to the 
measured value) and a B{E2)"^ of 0.163 e^b^ for the 7- 
vibrational state of ^^^Sm by solving the Schrodinger 
equation in collective quadrupole degrees of freedom 
(Bohr and Mottelson's collective model) with the Myers- 
Swiatechi potential. Garcia-Ramos et al. [68] used the in- 
teracting boson model (IBM) to obtain low-energy states 
in about 20 even-even rare-earth nuclei, eight of which 
are discussed here. For each isotopic chain they deter- 
mined the parameters of the IBM Hamiltonian by ap- 
proximately reproducing the measured excitation ener- 
gies for J- = 2+ 4+ 6+ 8+ 0+ 2+ 4+ 2+, 3+ and 
states, and determined the boson effective charge by 
reproducing several measured B{E2ys. With regard to 
the 7-vibrational B{E2)"\'^s of the eight nuclei computed 
here, they reproduced those of ^^^'^^^Gd well but overes- 
timated others. See also Ref. [69], which presented an- 
other set of IBM calculations. 



IV. /3-VIBRATIONS 

In Tab. Ill we show calculated and measured energies 
and B{E2)Ys, with SkM*, for "/^-vibrational" states. 
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FIG. 8: (Color online) Results of two calculations with different values for ^cut and the pairing strength Gq, where q=p (proton) 
or n (neutron). Calculation A (left) uses (^cut, Cn, Gp) = (60 MeV, 218.521 MeV fm^ 176.364 MeV fm^), and calculation B 
(right) uses (200 MeV, 168.384 MeV fm^, 143.996 MeV fm^). The SkM*-based results discussed in prior sections were obtained 
from Calculation A. 



As mentioned earlier, the = 0+ channel contains a 
spurious state, and we display only those "/3-vibrations" 
that are clearly uncontaminated by spurious motion (see 
the last column of the table). 



TABLE III: Properties of vibrational" states in four nuclei. 
-^cai ^fxp are the calculated (with SkM*) and experimen- 
tal energies. B(^2)tli and B(^2)tfxp are the corresponding 
reduced upward transition probabilities. The "Ratio of S'at" 
is the ratio of the (spurious) strength associated with the 
particle- number operator [10] for the vibrational" state to 
that for the spurious state (average of proton and neutron). 



Nucleus 




-'^exp 






Ratio of Sn 




(MeV) 


(MeV) 








166Yb 


1.802 


1.043 


0.0398 




0.004 


168Yb 


2.039 


1.155 


0.0343 




0.012 


172Yb 


1.605 


1.117 


0.0049 


0.0081(17) 


0.054 


i^°Er 


1.596 


0.960 


0.0030 


0.0079(9) 


0.054 



We obtained these results with calculation B (see 
above), the large cutoff in which should lead to a more 
accurate treatment of the spurious state, though contam- 
ination in nuclei not shown in the table indicates that a 
finer mesh is necessary with a large cutoff. In our previ- 
ous paper [10], which used ^cut = 60 MeV, E^^^ and 
5(^2) tf^i were 1.390 MeV and 1.117 e^b^ in ^^^Yb; 
our new energy is thus 15% larger. Compared to the 
7- vibrational states, overall we apparently overestimate 



energies and underestimate 5(£^2)t's, and do not do as 
good a job as with 7-vibrational sates. Reference [13] 
points out that "/^-vibrational" states are not purely vi- 
brational, and in many cases are better interpreted as 
the second member of the i^^=0+ yrare rotational band. 
The QRPA cannot describe rotational bands and so the 
discrepancy between our results and experiment is not 
totally surprising. 



V. CONCLUSION 

We have used the QRPA with the Skyrme functionals 
SkM* and SLy4 and volume-(5-pairing to calculate the 
energies and B{E2) t's of 7-vibrational states in well- 
deformed even-even rare-earth nuclei. SkM* proves to 
be the better functional. The range of calculated values 
overlaps well with that of the experimental data. Since 
the QRPA energies are appreciably different from their 
unperturbed counterparts, that counts as a success for 
the residual interaction; the vibrational states discussed 
here are not taken into account at all in determining en- 
ergy functionals. In detail, however, the calculations are 
far from perfect, and their TV- and Z- dependence sug- 
gest the importance of many-body correlations that are 
not included in the QRPA. And our representation of 
"/^-vibrational" states turns to be worse than that of 7 
vibrations, probably because "/3 vibrations" are often not 
really vibrations. 

We also suggested that the cutoff associated with vol- 
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ume pairing can be fixed along with the pairing strength 
by examining properties that are sensitive to the struc- 
ture of low-lying quasiparticles. 

Our calculation is better overall than the works of half 
a century ago. The aims of the pairing-plus-QQ model 
are much more limited than those of nuclear DFT; the 
mean field arising from pairing-plus-QQ is an infinitely 
deep well, and so the model cannot make predictions for 
binding energies or for excitation energies near the drip 
line (where the underlying Nilsson single-particle poten- 
tial is not appropriate). Despite the increasing sophisti- 
cation of the many-body methods that have been applied 
together with the pairng-plus-QQ model, a more general 
framework such as DFT appears necessary for the unified 
description of heavy nuclei. 

Finally, we have shown that in this era of supercom- 
puting a scalable code makes systematic and fully self- 
consistent Skyrme-QRPA studies possible. We expect, as 
a result, that excited states will play an increasing role 
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in the determination of nuclear density functionals. 
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Appendix: Coulomb-direct matrix elements 

The computation of the direct two-body matrix ele- 
ments of the Coulomb interaction consumes a lot of com- 
puting time. In this appendix we present our implemen- 
tation of that computation. 

The Coulomb interaction is 



\ri - r2\ 



(A.l) 



We take advantage of axial symmetry to write the wave 
function as 

X,(r) = -^ E •^a(a;^,p)e^(^«-<^)*|cT), (A.2) 
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The label a stands for (gTrj^i), i. e. particle type, par- 
ity, angular- momentum z-component, and an additional 
label i to fully specify the state. The position r is repre- 
sented in cylindrical coordinates, and the label a = ±1/2 
is the 2:-component of the spin. The function J^a(cr; p) 
is treated numerically. The set {Xa{r)} can refer to any 
single-particle basis (or components of quasiparticle basis 
states, in which case another label to distinguish upper 
from lower is necessary) with axial and parity symme- 
tries. In our calculations we use the canonical single- 
particle basis. 

With the help of a few well-known formulae from Ap- 
pendix B of Ref. [70] , one can obtain the expansion 



where 



pI + 4 
Pi + 4 



p% 
pI 



= pi + zi 

= pI + zI 



} , if p1 + zl <pI + zl , (A.4) 



1 



ri - r2 



E 



E 



{l-m)\ 



(A.3) 



and the Pim are associated Legendre polynomials [70]. 
By using Eqs. (A.l) — (A.4), one can then write the matrix 
element of the Coulomb-direct interaction as 



Valcd = J d'r^ J dV2Xt(ri)Xt(r2)Vc(ri,r2)Xe(ri)Xd(r2) 

^-i5-m+j-,0^-i-+m+i-,0 / dZi / dpipi / dZ2 / dp2p2J^a{cra]Zi,pi) 



(l-m)\ 



(A.5) 



(/-hm)! 
Changing variables to 

(z,i?= + , (A.6) 

and noting that 

Ta{(ya\ -Z, p) = (-)^'^"''^7ra ^^(cra; Z, p) , (A.7) 

Plm{-X) = {-y-^Plm{x) , (A.8) 

we arrive at 

Vab.cd = 4e^'^-J5+JJJf-jJ'^T„7re,7r67ri ^ ^ ^7r„7rc,(-)' L , / '^-^l^^Li^ / <^-R2-R2"^^72(-R2) 

+ /""di?iri(i?i)i?l+2 /"di?2^^^§^}| , (A.9) 



where 



Ti{Ri) = dziTa(aa; Zl, ^J Rj - zf^Pirn(^^^ , 

r2(i?2) = ^ ^""^ ^^2^6(^6; ^2, ^JR^ - ztjPim ^2, ^ Rl - zi) . (A.IO) 



Though it is not explicit in the notation, Ti{Ri) depends on the labels a, c, cTa, and {l,m), and 72(^^2) on similar 
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quantities. 

Equations (A. 9) and (A. 10) are what we use, with mi- 
nor modifications for hole states, in our code. We calcu- 
late Ti{Ri) and T2{R2) on a mesh and store them in ar- 
rays. Once this is finished, the time to calculate Eq. (A. 9) 
is determined mainly by the nest structure of the two-fold 
integrals and the summation with respect to /. For a sys- 
tem with quadrupole deformation y5 ~ 0.3 the number of 
terms necessary in the sum over / (much fewer than 20 
in practice, with only even or only odd / contributing) 
is much smaller than the number of mesh points in the 
integration. 



If an equidistant mesh is used for integrals in which 
an upper or lower bound is a variable, the computational 
effort to calculate the two-fold integrals in Eq. (A. 9) is 
nearly the same as that of single integrals. Thus, we 
calculate the wave functions on a new mesh by interpo- 
lating between B-spline points, and then use Simpson's 
rule with three times more mesh points than B-spline 
points to preserve accuracy (while still speeding up the 
integration). We have checked our procedure by using 
the two-body matrix elements that it produces to calcu- 
late the Coulomb-direct energy of the HFB ground state, 
which we then compared to the output of the HFB code. 



[1] S. C. Pieper, R. B. Wiringa, and J. Carlson, Phys. Rev. 

C 70, 054325 (2004). 
[2] S. C. Pieper and R. B. Wiringa, Annu. Rev. Nucl. Part. 

Sci. 51, 53 (2001). 
[3] A. G. Negoita, J. P. Vary, and S. Stoica, J. Phys. G 37, 

055109 (2010). 

[4] J. P. Vary, S. Popescu, S. Stoica, and P. Navratil, J. Phys. 

G 36, 085103 (2009). 
[5] G. Hagen, D. J. Dean, M. Hjorth- Jensen, T. Papenbrock, 

and A. Schwenk, Phys. Rev. C 76, 044305 (2007). 
[6] D. Dean and M. Hjorth-Jensen, Phys. Rev. C 69, 054320 

(2004). 

[7] J. Dobaczewski, e-print arXiv:nucl-th/1009.0899 (2010), 

to be published in J. of Phys. : Conference Series. 
[8] M. Kortelainen, T. Lesinski, J. More, W. Nazarewicz, 

J. Sarich, N. Schunck, M. V. Stoitsov, and S. Wild, Phys. 

Rev. C 82, 024313 (2010). 
[9] S. Goriely and E. Khan, Nucl. Phys. A 706, 217 (2002). 
[10] J. Terasaki and J. Engel, Phys. Rev. C 82, 034326 (2010). 
[11] S. Peru, G. Gosselin, M. Martini, M. Dupuis, S. Hilaire, 

and J.-C. Devaux, Phys. Rev. C 83, 014314 (2011). 
[12] K. Yoshida and T. Nakatsukasa, Phys. Rev. C 83, 

021304(R) (2011). 
[13] P. E. Garrett, J. Phys. G 27, Rl (2001). 
[14] M. Oshima, T. Morikawa, Y. Hatsukawa, S. Ichikawa, 

N. Shinohara, M. Matsuo, H. Kusakari, N. Kobayashi, 

M. Sugawara, and T. Inamura, Phys. Rev. C 52, 3492 

(1995). 

[15] M. Sugawara, H. Kusakari, T. Morikawa, H. In- 
oue, Y. Yoshizawa, A. Virtanen, M. Piiparinen, and 
T. Horiguchi, Nucl. Phys. A 557, 653 (1993). 

[16] C. Fahlander, B. Varnestig, A. Backlin, L. E. Svensson, 
D. Disdier, L. Kraus, I. Linck, N. Schulz, and J. Pedersen, 
Nucl. Phys. A 541, 157 (1992). 

[17] C. Fahlander, I. Thorslund, B. Varnestig, A. Backlin, 
L. E. Svensson, D. Disdier, L. Kraus, I. Linck, N. Schulz, 
J. Pedersen, et al., Nucl. Phys. A 537, 183 (1992). 

[18] B. Kotlinski, D. Cline, A. Backlin, K. G. Helmer, A. E. 
Kavka, W. J. Kernan, E. G. Vogt, C. Y. Wu, R. M. 
Diamond, A. O. Macchiavelli, et al., Nucl. Phys. A 517, 
365 (1990). 

[19] D. G. Burke, G. L0vh0iden, and T. F. Thorsteinsen, 

Nucl. Phys. A 483, 221 (1988). 
[20] T. Ichihara, H. Sakaguchi, M. Nakamura, M. Y. M. leiri, 

Y. Takeuchi, H. Togawa, T. Tsutsumi, and S. Kobayashi, 

Phys. Rev. C 36, 1754 (1987). 



[21] P. M. Walker, Phys. Scr. T5, 29 (1983). 

[22] R. M. Ronningen, R. S. Grantham, J. H. Hamilton, R. B. 

Piercey, A. V. Ramayya, B. van Nooijen, H. Kawakami, 

W. Lourens, R. S. Lee, W. K. Dagenhart, et al., Phys. 

Rev. C 26, 97 (1982). 
[23] J. R. Cresswell, P. D. Forsyth, D. G. E. Martin, and R. C. 

Morgan, J. Phys. G 7, 235 (1981). 
[24] F. K. McGowan and W. T. Milner, Phys. Rev. C 23, 

1926 (1981). 

[25] L. L. Riedinger, E. G. Funk, J. W. Mihelich, G. S. 

Schilling, A. E. Rainis, and R. N. Oehlberg, Phys. Rev. 

C 20, 2170 (1979). 
[26] F. K. McGowan, W. T. Milner, R. L. Robinson, P. H. 

Stelson, and Z. W. Grabowski, Nucl. Phys. A 297, 51 

(1978). 

[27] H. J. Wollersheim and Th. W. Elze, Z. Phys. A 280, 277 
(1977). 

[28] R. M. Ronningen, J. H. Hamilton, A. V. Ramayya, 

L. Varnell, G. Garcia-Bermudez, J. Lange, W. Lourens, 

L. L. Riedinger, R. L. Robinson, P. H. Stelson, et al., 

Phys. Rev. C 15, 1671 (1977). 
[29] C. W. Reich, R. C. Greenwood, and R. A. Lokken, Nucl. 

Phys. A 228, 365 (1974). 
[30] C. Baktash, J. X. Saladin, J. O'Brien, I. Y. Lee, and J. E. 

Holden, Phys. Rev. C 10, 2265 (1974). 
[31] R. N. Oehlberg, L. L. Riedinger, A. E. Rainis, A. G. 

Schmidt, E. G. Funk, and J. W. Mihelich, Nucl. Phys. A 

219, 543 (1974). 
[32] M. H. Cardoso, P. F. A. Goudsmit, and J. Konijn, Nucl. 

Phys. A 205, 121 (1973). 
[33] C. E. Bemis Jr., P. H. Stelson, F. K. McGowan, W. T. 

Milner, J. L. C. Ford Jr., R. L. Robinson, and W. Tuttle, 

Phys. Rev. C 8, 1934 (1973). 
[34] J. M. Domingos, G. D. Symons, and A. C. Douglas, Nucl. 

Phys. A 180, 600 (1972). 
[35] M. T. Gillin and N. F. Peek, Phys. Rev. 4, 1334 (1971). 
[36] A. Charvet, D. H. Phuoc, R. Duffait, A. Emsallem, and 

R. Chery, J. de Physique 32, 359 (1971). 
[37] H. Ejiri and G. B. Hageman, Nucl. Phys. A 161, 449 

(1971). 

[38] T. Grotdal, K. Nyb0, T. Thorsteinsen, and B. Elbek, 

Nucl. Phys. A 110, 385 (1968). 
[39] E. Veje, B. Elbek, B. Herskind, and M. C. Olesen, Nucl. 

Phys. A 109, 489 (1968). 
[40] R. Bloch, B. Elbek, and P. O. Tj0m, Nucl. Phys. A 91, 

576 (1967). 



11 



[41] Y. Yoshizawa, B. Elbek, B. Herskind, and M. C. Olesen, 

NucL Phys. 73, 273 (1965). [56 
[42] O. Nathan and V. L Popov, NucL Phys. 21, 631 (1960). [57 
[43] [http: / / www.nndc.bnl.gov] . 

[44] J. Bartel, P. Quentin, M. Brack, C. Guet, and H.-B. 

Hakansson, Nucl. Phys. A 386, 79 (1982). [58 
[45] E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and 

R. Schaeffer, NucL Phys. A 635, 231 (1998). [59 
[46] J. Terasaki, J. Engel, M. Bender, J. Dobaczewski, 

W. Nazarewicz, and M. Stoitsov, Phys. Rev. C 71, [60" 

034310 (2005). 

[47] M. V. Stoitsov, J. Dobaczewski, W. Nazarewicz, S. Pittel, [61 
and D. J. Dean, Phys. Rev. C 68, 054312 (2003). [62 

[48] A. Blazkiewicz, V. E. Oberacker, A. S. Umar, and 

M. Stoitsov, Phys. Rev. C 71, 054321 (2005). [63 

[49] CD. Boor, A Practical Guide to Splines (Springer, New 

York, 1978). [64 

[50] G. Niirnberger, Approximation by Spline Functions 

(Springer, New York, 1989). [65 

[51] L. L. Schumaker, Spline Function : Basic Theory (Cam- 
bridge Univ. Press, Cambrige, 2007). [66 

[52] A. Bohr and B. R. Mottelson, Nuclear Structure, voL 1 

(Benjamin, New York, 1969). [67 

[53] J. Terasaki, J. Engel, and G. F. Bertsch, Phys. Rev. C [68 
78, 044311 (2008). 

[54] J.-P. Delaroche, M. Girod, J. Libert, H. Goutte, S. Hi- [69 
laire, S. Peru, N. Pillet, and G. F. Bertsch, Phys. Rev. C 
81, 014303 (2010). [70 

[55] G. F. Bertsch, M. Girod, S. Hilaire, J.-P. Delaroche, 
H. Goutte, and S. Peru, Phys. Rev. Lett. 99, 032502 



(2007). 

A. Bulgac and Y. Yu, Phys. Rev. Lett. 88, 042504 (2002). 
V. O. Nesterenko, W. Kleinig, J. Kvasil, P. Vesely, P.- 
G. Reinhard, and D. S. Dolci, Phys. Rev. C 74, 064306 
(2006). 

A. P. Severyukhin, V. V. Voronov, and N. V. Giai, Phys. 
Rev. C 77, 024322 (2008). 

E. R. Marshalek and J. O. Rasmussen, Nucl. Phys. 43, 
438 (1963). 

D. R. Bes, P. Federman, E. Maqueda, and A. Zuker, Nucl. 
Phys. 65, 1 (1965). 

M. Zielinska-Pfabe, Act. Phys. Pol. B 2, 207 (1971). 

L Hamamoto, Suppl. Prog. Theor. Phys. 74 and 75, 157 

(1983). 

M. Matsuo and K. Matsuyanagi, Prog. Theor. Phys 78, 
591 (1987). 

T. Kishimoto and T. Tamura, Nucl. Phys. A 270, 317 
(1976). 

V. G. Soloviev and N. Y. Shirikova, Z. Phys. A 334, 149 
(1989). 

V. G. Soloviev and A. V. Sushkov, Z. Phys. A 345, 155 
(1993). 

K. Kumar, Nucl. Phys. A 92, 653 (1967). 

J. E. Garcia-Ramos, J. M. Arias, J. Barea, and A. Frank, 

Phys. Rev. C 68, 024307 (2003). 

D. D. Warner and R. F. Casten, Phys. Rev. C 28, 1798 
(1983). 

A. Messiah, Quantum Mechanics (North Holland, Ams- 
terdam, 1961). 



